#include "GridGenerator.h"
#include <fstream>
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <limits>

GridGenerator::GridGenerator(const SurfaceExtractor& surfaces, double size)
    : surfaceExtractor(surfaces), gridSize(size) {
    logFile.open("grid_generator_log.txt");
    if (!logFile) {
        throw std::runtime_error("无法创建网格生成器日志文件");
    }
}

void GridGenerator::computeSpaceRange() {
    // 获取所有顶点
    const auto& vertices = surfaceExtractor.getVertices();
    
    if (vertices.empty()) {
        throw std::runtime_error("没有找到任何顶点数据");
    }

    // 初始化范围为第一个点的坐标
    range.x_min = range.x_max = vertices[0].x;
    range.y_min = range.y_max = vertices[0].y;
    range.z_min = range.z_max = vertices[0].z;

    // 遍历所有顶点找到范围
    for (const auto& vertex : vertices) {
        range.x_min = std::min(range.x_min, vertex.x);
        range.x_max = std::max(range.x_max, vertex.x);
        range.y_min = std::min(range.y_min, vertex.y);
        range.y_max = std::max(range.y_max, vertex.y);
        range.z_min = std::min(range.z_min, vertex.z);
        range.z_max = std::max(range.z_max, vertex.z);
    }

    logFile << "空间范围：" << std::endl;
    logFile << "X: [" << range.x_min << ", " << range.x_max << "]" << std::endl;
    logFile << "Y: [" << range.y_min << ", " << range.y_max << "]" << std::endl;
    logFile << "Z: [" << range.z_min << ", " << range.z_max << "]" << std::endl;
}

void GridGenerator::generateNodes() {
    // 计算每个方向的网格数量和实际网格尺寸
    double x_length = range.x_max - range.x_min;
    double y_length = range.y_max - range.y_min;
    double z_length = range.z_max - range.z_min;

    // 计算每个方向的网格数量，向下取整以确保实际网格尺寸大于基础尺寸
    int nx = static_cast<int>(std::floor(x_length / gridSize));
    int ny = static_cast<int>(std::floor(y_length / gridSize));
    int nz = static_cast<int>(std::floor(z_length / gridSize));

    // 确保至少有一个网格
    nx = std::max(1, nx);
    ny = std::max(1, ny);
    nz = std::max(1, nz);

    // 计算实际的网格尺寸（会大于等于基础网格尺寸）
    actual_dx_ = x_length / nx;
    actual_dy_ = y_length / ny;
    actual_dz_ = z_length / nz;

    logFile << "网格分辨率：" << nx << " x " << ny << " x " << nz << std::endl;
    logFile << "实际网格尺寸：" << std::endl;
    logFile << "dx = " << actual_dx_ << std::endl;
    logFile << "dy = " << actual_dy_ << std::endl;
    logFile << "dz = " << actual_dz_ << std::endl;

    // 生成规则网格点
    nodes.clear();
    nodes.reserve((nx + 1) * (ny + 1) * (nz + 1));

    for (int k = 0; k <= nz; ++k) {
        for (int j = 0; j <= ny; ++j) {
            for (int i = 0; i <= nx; ++i) {
                Point3D node;
                node.x = range.x_min + i * actual_dx_;
                node.y = range.y_min + j * actual_dy_;
                node.z = range.z_min + k * actual_dz_;
                node.type = 0;  // 普通节点
                nodes.push_back(node);
            }
        }
    }

    logFile << "生成节点数量：" << nodes.size() << std::endl;
}

void GridGenerator::generateTetrahedra() {
    // 使用相同的网格数量计算
    double x_length = range.x_max - range.x_min;
    double y_length = range.y_max - range.y_min;
    double z_length = range.z_max - range.z_min;

    int nx = static_cast<int>(std::floor(x_length / gridSize));
    int ny = static_cast<int>(std::floor(y_length / gridSize));
    int nz = static_cast<int>(std::floor(z_length / gridSize));

    nx = std::max(1, nx);
    ny = std::max(1, ny);
    nz = std::max(1, nz);

    elements.clear();
    elements.reserve(nx * ny * nz * 5);  // 每个立方体分成5个四面体

    // 生成四面体单元
    for (int k = 0; k < nz; ++k) {
        for (int j = 0; j < ny; ++j) {
            for (int i = 0; i < nx; ++i) {
                // 计算立方体的8个顶点索引
                int nx1 = nx + 1;
                int ny1 = ny + 1;
                size_t v000 = k * nx1 * ny1 + j * nx1 + i;
                size_t v100 = v000 + 1;
                size_t v010 = v000 + nx1;
                size_t v110 = v010 + 1;
                size_t v001 = v000 + nx1 * ny1;
                size_t v101 = v001 + 1;
                size_t v011 = v001 + nx1;
                size_t v111 = v011 + 1;

                // 使用Kuhn剖分方案，确保四面体的正确方向
                // 每个四面体的顶点顺序遵循右手法则
                elements.push_back({v000, v110, v011, v010});  // 底面逆时针，顶点在上
                elements.push_back({v000, v101, v110, v100});
                elements.push_back({v000, v011, v101, v001});
                elements.push_back({v110, v101, v011, v111});
                elements.push_back({v000, v110, v101, v011});
            }
        }
    }

    logFile << "生成四面体单元数量：" << elements.size() << std::endl;
}

bool GridGenerator::generateInitialGrid() {
    try {
        logFile << "开始生成初始网格..." << std::endl;
        logFile << "基础网格尺寸：" << gridSize << std::endl;

        // 计算空间范围
        computeSpaceRange();

        // 生成节点
        generateNodes();

        // 生成四面体单元
        generateTetrahedra();

        logFile << "初始网格生成完成" << std::endl;
        return true;
    }
    catch (const std::exception& e) {
        logFile << "错误：" << e.what() << std::endl;
        return false;
    }
}

void GridGenerator::writeToVTK(const std::string& filename) {
    std::ofstream file(filename);
    if (!file) {
        throw std::runtime_error("无法创建VTK文件：" + filename);
    }

    // 写入VTK文件头
    file << "# vtk DataFile Version 3.0\n";
    file << "Initial Tetrahedral Grid\n";
    file << "ASCII\n";
    file << "DATASET UNSTRUCTURED_GRID\n\n";

    // 写入节点
    file << "POINTS " << nodes.size() << " double\n";
    for (const auto& node : nodes) {
        file << std::fixed << std::setprecision(6)
             << node.x << " " << node.y << " " << node.z << "\n";
    }
    file << "\n";

    // 写入单元
    file << "CELLS " << elements.size() << " " << elements.size() * 5 << "\n";
    for (const auto& element : elements) {
        file << "4 " << element.v1 << " " << element.v2 << " "
             << element.v3 << " " << element.v4 << "\n";
    }
    file << "\n";

    // 写入单元类型
    file << "CELL_TYPES " << elements.size() << "\n";
    for (size_t i = 0; i < elements.size(); ++i) {
        file << "10\n";  // VTK_TETRA = 10
    }
    file << "\n";

    file.close();
    logFile << "VTK文件已保存：" << filename << std::endl;
} 